Statistical multi-moment bifurcations in random delay coupled swarms 
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We study the effects of discrete, randomly distributed time delays on the dynamics of a coupled 
system of self-propelling particles. Bifurcation analysis on a mean field approximation of the system 
reveals that the system possesses patterns with certain universal characteristics that depend on 
distinguished moments of the time delay distribution. Specifically, we show both theoretically and 
numerically that although bifurcations of simple patterns, such as translations, change stability only 
as a function of the first moment of the time delay distribution, more complex patterns arising from 
Hopf bifurcations depend on all of the moments. 



Recently, much attention has been given to the study 
of interacting mult i- agent, particle or swarming systems 
in various natural and engineering fields. Interestingly, 
these multi-agent swarms can self-organize and form 
complex spatio-temporal patterns even when the cou- 
pling between agents is weak. Many of these investi- 
gations have been motivated by a multitude of biological 
systems such as schooling fish, swarming locusts, flock- 
ing birds, bacterial colonies, ant movement, etc. p]-|U, 
and have also been applied to the design of systems of 
autonomous, communicating robots or agents 0-01 and 
mobile sensor networks ||. 

Many studies describe the swarm system at the in- 
dividual, or particle, level via models constructed with 
ordinary differential equations (ODEs) or delay differen- 
tial equations (DDEs) to describe the trajectories j9l[Io|. 
When there are a large number of densely-distributed 
particles, authors have employed partial differential equa- 
tions (PDEs) to describe the average agent density and 
velocity 0, [j, [HI, [H[. Recently, the inclusion of noise 
in such particle-based studies has revealed interesting, 
noise-induced transitions between different coherent pat- 
terns [H, [3 ■ Such noise driven systems have led to the 
discovery of first and second order phase transitions in 
swarm models [l5j . 

A topic of intense ongoing research in interacting parti- 
cle systems, and in particular in the dynamics of swarms, 
is the effect of time delays. It is well known that time 
delays can have profound dynamical consequences, such 
as destabilization and synchronization [lH Il7|. and de- 
lays have been effectively used for purposes of control 
fl8| . Initially, such studies focused on the case of one or 
a few discrete time delays. More recently, however, the 
complex situation of several and random time delays has 
been researched [l9l-|2lj]. An additional important case 
is that of distributed time delays, when the dynamics of 
the system depends on a continuous interval in its past 
instead of on a discrete instant [12, HH . 

There exists a complex interplay between the attrac- 
tive coupling, time delay, and noise intensity that pro- 
duces transitions between different spatio-temporal pat- 
terns [13, HH in the case of a single, discrete delay. 
Here, we consider a more general swarming model where 
coupling information between particles occurs with ran- 



domly distributed time delays. We perform a bifurcation 
analysis of a mean field approximation and reveal the pat- 
terns that are possible at different values of the coupling 
strength and parameters of the time delay distribution. 

We model the dynamics of a 2D system of N identical 
self-propelling agents that are attracted to each other in 
a symmetric manner. We consider the effects of finite 
communication speeds and information-processing times 
so that the attraction between agents occurs in a time 
delayed fashion. The time delays are non-uniform but 
they are symmetric among agents Ty(= Tji), for parti- 
cles i and j, as well as constant in time. The dynamics 
of the particles is described by the following governing 
equations: 
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for i = 1,2 . . . , N. The vector r$ denotes the position of 
the ith agent at time t. The term (l — |i"i| 2 ) i"i represents 
self-propulsion and frictional drag forces that act on each 
agent. The coupling constant a measures the strength of 
the attraction between agents and the time delay between 
particles i and j is given by r,j . When a = the agents 
tend to move in a straight line with unit speed as time 
tends to infinity. The N(N — l)/2 different time delays 
< t,j(= Tji) are drawn from a distribution p(r) whose 
mean and standard deviation are denoted by fx T and a T , 
respectively. 

We obtain a mean field approximation of the swarming 
system by measuring the particle's coordinates relative 
to the center of mass = R + <5r,, for i = 1,2 ... ,N, 
where R(f) = ^X)ili r i(0- Following the approxima- 
tions from [25l ] , we obtain a mean field description of the 
swarm: 
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The approximations 
© require that N 



necessary to obtain Eq. 
be sufficiently large so that 

) « / °°R(t-r)p(r)dr 



N(N-l) s—ii=L ^} 

and that the swarm particles remain close together 
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FIG. 1. Bifurcation structure of the translating state of the 
mean field Eqs. (J2J) for the exponential distribution described 
in the text (a T = 0.2). The translating state merges with: (i) 
the stationary state along the continuous red curve ap T = 1; 
and (ii) the circularly rotating state along the dashed black 
curve a(r 2 } = 2. The component of the translating state 
parallel to the motion undergoes a Hopf bifurcation along the 
green, dotted curve. (Color online.) 



Equations ^ admit a uniformly translating solution 
R(i) = Ro + Vo • t (Ro and Vo are any constant 2D 
vectors). The speed |V | must satisfy 
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which shows that this solution is possible as long as the 
system parameters lie below the hyperbola ap T = 1 in 
the (a, p T ) plane. Remarkably, the speed of the of the 
translating state depends exclusively on the mean of the 
distribution p{r) and not on any of the higher moments. 

The linear stability of the translating state is examined 
by taking X(t) = y/1 - ap T -t + 5X(t) and Y(t) = SY(t). 
The two linearized equations decouple and the stability 
of motions parallel and perpendicular to the translating 
direction are determined by the characteristic equations 
Dm (A) and Dj_(A), respectively: 



Dm (A) = J"(A) - (3ap T - 2) A, V± (A) = T{\) 



ap T X, 
(4) 
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where J"(A) = a (l - (e~ Ar )) + A 2 . The function 
is the moment generating function of p(r) since the n-th 
moment is equal to (r n ) = (-l)"^r(e- AT )| A= o- Re- 
gardless of the choice of a and p(r), the characteristic 
functions Dy and V± have a zero eigenvalue arising from 
the translation invariance of Eq. ((2J H(| . There is a fold 
bifurcation as an eigenvalue of Dm crosses the origin when 
ap T = 1 , which marks the disappearance of the translat- 
ing state as seen from Eq. ©. Numerical analysis [27j 
reveals an additional curve on the (a, p T ) plane (below 
the curve ap T = 1) along which perturbations parallel to 
the translation undergo a Hopf bifurcation as a complex 
pair of eigenvalues of Dp cross the imaginary axis. 

As for perturbations perpendicular to the translational 
motion, there is another fold bifurcation as an eigenvalue 
of V± crosses the origin along the curve a(r 2 ) = 2, which 
represents a bifurcation in which the translating state 
merges with a circularly rotating state of infinite radius, 
as discussed below. 

Considering a fixed er T , the overall stability picture of 
the translating state of Eqs. is as follows (see Fig. 



[T]). For values of (a, p T ) below the curves a(r) = 1 and 
a(r 2 ) = 2 (region A) the translating state is linearly sta- 
ble . These two curves may cross at a point that we 
call the 'zero frequency Hopf point' (ZFH). The trans- 
verse direction of the translating state becomes unstable 
along the curve a(r 2 ) = 2 where this state merges with 
the circularly rotating state (along the mentioned curve 
the rotating state has an infinite radius); transverse per- 
turbations of the translating state will thus produce a 
transition to the rotating state in regions B and C. From 
the ZFH point, there emanates a Hopf bifurcation curve 
where the parallel component of the translating state be- 
comes unstable so that in region C there is a transition 
from the translating state to oscillations along a straight 
line. Finally, the translating state ceases to exist along 
the curve o(r) = 1 where there is a pitchfork type bi- 
furcation with the stationary steady state solution. The 
possible behaviors in region D are discussed below. 
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FIG. 2. Speed of the translating state of Eqs. |T} and ((2)1 as 
a function of p T ; here N = 150, a = 0.2 (a) and a = 1.2 (b). 
The red line (color online) represents the mean field result 
from Eq. (J2|; the continuous segment marks where the trans- 
lating state is linearly stable and the dashed segment where 
it is unstable. The symbols represent numerical simulations 
of Eqs. fT]) for different time delay distributions: sliding ex- 
ponential and sliding uniform o> — 0.5 (a), a T = 0.05 (b); 
widening exponential p T — <j t (a) and (b); widening uniform 
Pt = V3(T r (a) and (b). 

We compare the mean field bifurcation results with the 
full swarm system via numerical simulations. Here, we 
make use of two different time delay distributions with 
mean p T and standard deviation <r r to test our find- 
ings. The first is an exponential distribution p(r) = 

e /<7 T for t > p T — a T and zero otherwise; we 

require a T < p T for proper normalization. The second 
distribution is a uniform p(r) = ^= - for p T — \/Za T < 

t < p T + v3er T and zero otherwise; here, we require 
< p T . Moreover, we employ two versions of the 
mentioned distributions: a 'sliding' one (er T = const.) 
and a 'widening' version (/i r cx a T ). 

Fig. [2] compares the speed of the swarm obtained from 
Eq. §5§ with the time-averaged speed of the center of 
mass obtained from simulations (after the decay of tran- 
sients). The swarm particles arc all located at the origin 
at time zero and move with the speed obtained from Eq. 
(O along the x axis. In these simulations, we use both 
'sliding' (a> fixed) and 'widening' (both p T and o T vary) 
versions of the exponential and uniform distributions de- 
scribed above. Fig. [2] shows that the swarm converges to 
the translating state up to a value of p T beyond which 
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the swarm converges to a state in which it oscillates back 
and forth along a line with a near zero time-averaged 
speed (the average is taken over an interval much longer 
than the period). The transition to the oscillatory regime 
occurs earlier than the mean field prediction. The full 
simulation results show that in this deviation from the 
mean field, the swarm particles become spread out too 
far apart and render the approximations leading to Eq. 
© invalid. 

10 



\ c \ '• \ 


a 


\ \ B '• ' • 









10 



10 



10 





FIG. 3. Bifurcation curves of the mean field Eqs. ([2| at 
fixed ay for the two time delay distributions p(r) described 
in the text: exponential (top) and uniform (bottom). The 
translational state disappearance curve ap T — I (red), bifur- 
cation of the translational state with circularly rotating state 
curve o(t 2 } = 2 (black). The first four members of the sta- 
tionary state Hopf bifurcation curves are also shown (blue, 
dashed green, dotted-dashed cyan and dotted magenta). In 
each panel, ay has the value (a) 0.2, (b) 0.95, (c) 0.2 and (d) 
0.3, respectively. (Color online.) 

In addition to the translating state, Eqs. ([2]) always 
possess a stationary state solution R(£) = R = const. 
In the full system, Eq. (UJ, the stationary state for the 
center of mass manifests itself in a swarm 'ring state', 
where some particles rotate clockwise and others counter- 
clockwise on a circle around a static center of mass. The 
characteristic equation that governs the linear stability 
of the stationary state has the form (25(A)) = 0, where 
T>(\) = J 7 (A) — A. Once more there is a zero eigenvalue 
for all choices of a and p(r) that arises from the trans- 
lation invariance of Eqs. ©. Also, since 2?(0) = 0, 
2?'(0) = a/i r — 1 and lim\-> 00 T>(X) = oo, the condition 
a/i T — 1 < guarantees the existence of at least one 
real and positive eigenvalue which renders the station- 
ary state linearly unstable. Thus, afJ, T = 1 is a bifurca- 
tion curve on the (a, fj, T ) plane along which the uniformly 
translating state bifurcates with the stationary state. 

The stationary state undergoes Hopf bifurcations when 
the equation V(iuj) = a (l — (e~'" T )) — uj 2 — ioj — for 
uj 7^ is satisfied. The function {e~ %U3T ) is called the char- 
acteristic function of p(r) and is related to the moment 
generating function of the distribution; its Taylor series 
contains all of the moments of p(r). This shows that the 



location of the Hopf bifurcations depends on the values 
of all moments of the time delay distribution. This is in 
contrast to the region where the translating state exists 
ap T < 1, which involves the first moment only. 

At the location of the Hopf bifurcations, circular orbits 
bifurcate from the stationary state. This may be seen by 
changing Eq. © from the Cartesian (X, Y) to polar 
coordinates (R, (f>), noticing that circular orbits R = Rq, 
<p = Lot are possible as long as 



a (1 — (coswr)) , 



Rq 



1 (sincjr), 



(5) 



and then realizing that the first of (JS|) and the condition 
Rq = are precisely the real and imaginary parts of the 
Hopf bifurcation conditions T>(iuj) = 0. 





FIG. 4. Center of mass speed as a function of fi T for the 'slid- 
ing' exponential time delay distribution (left) and the 'sliding' 
uniform distribution (right). Here N = 150, a = 2. The con- 
tinuous red curve represents the mean field result from Eq. 
((SJ), while the symbols represent the results from numerical 
simulations of Eqs. Q. (color online) 

Generically, the Hopf conditions for the stationary 
state T>(iuS) = yield a family of curves in the (a, p T ) 
plane (Fig. The first member of the Hopf family 

emanates from the crossing of the curves ap T = 1 and 
a(r 2 ) = 2; the former curve is where the translating state 
bifurcates from the stationary state in a pitchfork-like bi- 
furcation. Hence the name 'Zero Frequency Hopf for the 
Hopf- fold point (Fig. (3^)- The first Hopf curve is super- 
critical and gives rise to a circularly rotating orbit with 
radius and frequency given by the first solution of Eqs. 
©. Below this first Hopf curve and ap T — 1, in region 
A the stationary state is stable. From Eqs. ([5]) it follows 
that this circularly rotating orbit collides with the trans- 
lating state along the curve a(r 2 ) = 2 where its radius 
tends to infinity and its speed to that of the translating 
state, y/l — a\i T . Thus, in regions B and C the system 
converges to the circularly rotating orbit. The different 
regions change shape for the other panels of Fig. [3j but 
the dynamics remain as described above. 

We compare the mean field prediction for the location 
of the birth of the circularly rotating state (first Hopf 
curve) with the full system. Fig. |4] shows the results 
from numerical simulations of Eqs. ([T]) at a fixed value 
of ay for increasing fx T . We plot the speed of the cen- 
ter of mass averaged over a long time interval, after the 
decay of transients. In that plot, the near-zero values of 
the mean speed (in the interval < fi T < 1.4) indicate 
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that the particles have converged to the ring state, while 
for all higher values of p T the swarm converges to the 
rotating state. Remarkably, for the highest values of p T , 
the center of mass of the swarm moves faster than unit 
speed, the asymptotic speed of uncoupled particles. The 
reason is that while in its rotating orbit, the mean time- 
delayed position of the swarm is actually "ahead" of the 
center of mass at the current time, causing the particles 
to accelerate forward along the circular orbit. 

In summary, we have considered a randomly delay dis- 
tributed coupled swarm model, and analyzed the bifur- 
cations of various patterns as a function of delay char- 
acteristics and coupling strength. In particular, we have 
shown that the location and shape of the Hopf bifurca- 
tion curves is strongly-dependent on all the moments of 
p(r). This dependence, in addition to the fact that the 



succeeding Hopf curves in Fig. [3] exhibit higher frequen- 
cies of rotation, makes the higher-order patterns equally 
sensitive to all the moments of the delay distribution. 
In the single delay case with distribution p = <5(r — To), 
where S(t) is a Dirac delta function, all of the succeeding 
Hopf bifurcations are all subcritical and continuous. In 
contrast, when all moments are present, the bifurcations 
may not even be continuous, presenting their structure 
as isolated closed curves bounded by fold bifurcations, as 
seen in Fig.[3]i for a uniform distribution. Finally, we ex- 
pect that in other globally delay-coupled systems |28l - [31| , 
generic types of behavior involving bifurcations including 
all moments of the distribution should be present. 
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